Friedel oscillations at the Dirac-cone-merging point in anisotropic graphene 
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We study the Friedel oscillations induced by a localized impurity in anisotropic graphene. We 
focus on the limit when the two inequivalent Dirac points merge. We find that in this limit the 
Friedel oscillations manifest very peculiar features, such as a strong asymmetry and an atypical 
inverse square-root decay. Our calculations are performed using both a T-matrix approximation 
and a tight-binding exact diagonalization technique. They allow us to obtain numerically the local 
density of states as a function of energy and position, as well as an analytical form of the Friedel 
oscillations in the continuum limit. The two techniques yield results that are in excellent agreement, 
confirming the accuracy of such methods to approach this problem. 



I. INTRODUCTION 

Graphene has known an increased interest over the 
past years, with some of the most interesting questions 
at present focusing on the possibility to modify the elec- 
tronic structure of graphene, either by mechanical defor- 
mations such as stretching, or twistin^^, via chemical 
additions, or by changing the nature of substrate. In the 
light of possible important applications, the most promis- 
ing directions have been towards opening a gapP, enhanc- 
ing the spin-orbit interaction, the realization of the quan- 
tum spin Hall effect^ ^, and obtaining integer and frac- 
tional quantum Hall states using pseudo-magnetic (cur- 
vature) field^I^. 

One of the most studied modification of graphene is 
mechanical stretching, which gives rise to a hopping 
anisotropy, and consequently to a strong renormaliza- 
tion of the band structure. For a critical value of the 
anisotropy, such normalization has been predicted to give 
rise to a hybrid Dirac cone, exhibiting a linear disper- 
sion along one direction, and a quadratic one along the 
perpendicular one^. Recently, the realization of a cold- 
atom equivalent of such anisotropic system has been 
achieved^. Similar hybrid Dirac cones have been pre- 
dicted to arise when the higher-order hopping parameters 
are strongly enhanced^^^^, which may occur for example 
in the presence of adatom^. 

In this work we focus on a system with such hybrid 
semi-Dirac points and we study the Friedel oscillations 
(FO) generated in the presence of a single localized im- 
purity. We use both analytical techniques such as the 
T-matrix approximation, and numerical techniques (the 
exact diagonalization of the lattice tight-binding Hamil- 
tonian). Using the T-matrix approximation we obtain 
the form of the Fourier transform of the Friedel oscil- 
lations induced by the impurity. We also calculate the 
real-space form of these oscillations. For small energies 
and large distances (in the continuum limit) we obtain 
an exact analytical form of these oscillations, while we 
evaluate the short distance behavior of the Friedel oscil- 
lations using a numerical integration. On the other hand. 



we calculate the local density of states (LDOS) at each 
lattice site using an exact diagonalization of the lattice 
tight-binding Hamiltonian. Finally we study the form 
of the LDOS at zero energy using wavefunction argu- 
ments along the lines of Ref. [TH which allow us to obtain 
an analytical form for the impurity state at zero energy. 
The results obtained via the above methods are in per- 
fect agreement, confirming the accuracy of these tools for 
describing the impurity effects in such systems. 

The most interesting characteristic of the observed 
Friedel oscillations is a strong anisotropic spatial depen- 
dence - the period and decay length of these oscilla- 
tions depends strongly on direction - consistent with the 
anisotropy of the band structure. Also we observe an 
atypical inverse square root decay for large distances and 
small energies on each of the two A and B sublattices. 
Moreover, similar to the isotropic graphene, the LDOS 
contributions of the two sublattices are dephased by tt, 
yielding a cancelation of the 1 / ^/r terms and an effective 
1/r decay of these oscillations with the distance from the 
impurity. 

The structure of the paper is as follows. In section II we 
present the model employed to describe isotropic as well 
as anisotropic graphene. In Section III we present the 
Friedel oscillations in the LDOS calculated using wave- 
function considerations (A), tight-binding exact diago- 
nalization (B) and the T-matrix approximation (C). We 
conclude in section IV. 



II. MODEL 

Graphene consists of a honeycomb lattice of carbon 
atoms with two atoms {A and B) per unit cell (see Fig.jl]). 
Denoting the distance between two nearest neighbors 
ao, with ao = 0.142 nm, then ai= ao(— ^,|) and 

^2= ao(^, |) are basis vectors of the triangular Bra- 
vais lattice. 

The corresponding first Brillouin Zone (BZ) is hexag- 
onal, as depicted by the green dashed line in Fig. [2] Its 
geometrical properties only depend on the Bravais lat- 



FIG. 1. Graphene honeycomb lattice. 



tice. Nevertheless the number of atoms per unit cell be- 
comes relevant for the energy spectrum. For graphene 
(two atoms per unit cell with one electron per atom) the 
energy bands are well described using a tight-binding 
model: each 2pz electron may hop between two sites i 
and j with a given amplitude tij. In this work we only 
consider the hopping between nearest neighbors, with a 
fixed hopping amplitude t ~ 2.7 eV for the nearest neigh- 
bor vectors Ji, ^2, and a variable amplitude t' for 5s. The 
corresponding second-quantized Hamiltonian is given by: 
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and f{k) = -t{e-'^-^^ + ^-'^-32^ _ fe-'^-^K Sbz is the 
area of the BZ. The operators a {b) and (6'*') are field 
operators that respectively annihilate and create an elec- 
tron on the A (B) sublattice. The energy spectrum is 
then obtained by diagonalizing the Hamiltonian matrix 
As there are two atoms per unit cell, there are two 

energy bands e±{k) = Negative values of e cor- 

respond to the valence band whereas positive ones corre- 
spond to the conduction band. When t' = t there are two 
inequivalent points K and K' at the corners of the BZ 
for which the two bands touch. These points are denoted 
Dirac points since the energy spectrum is conical in their 
vicinity. Note that the coincidence between the Dirac 
points (determined by the band structure) and that of 
the corners of the BZ (intrinsic to the Bravais lattice) 
occurs only when f = t. 

Fig. |2] illustrates the fact that the Dirac points move 
away from the corners of the BZ when varying the hop- 
ping parameter t' . Increasing this amplitude fr om t to 
2t makes the two inequivalent Dirac points merge'^^"'^ at 
the M point (right at the middle of edges of the BZ). 
The critical value t' = 2t corresponds to the annihila- 
tion of a pair of Dirac points with opposite Berry phases. 
This topological invariant changes abruptly from ±7r to 
at the merging, which thus defines a topological transi- 
tion between a semi-metallic phase and a band insulator. 
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FIG. 2. Energy spectra when t' — t (top) and t' — 2t (bot- 
tom). The later corresponds to the merging of Dirac points 
into a single point M. The green dashed line depicts the 
Brillouin zone. 



since a gap opens at the M point for t' > 2t. This can 
be seen by expanding f{k) in the vicinity of this point 



defined by (0, |^): 



Here 
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and A = t' — 2t charac- 



terizes the distance from the topological transition and 
also gives the value of the gap when t' > 2t. Exactly at 
the transition (A = 0), the Hamiltonian exhibits a semi- 
Dirac energy dispersion such that (q) = ±\f^{q)\ is 
linear in qy but quadratic with respect to qx'. 



e^{q) = ±d{cyqyy 



2m* 



(4) 



III. LOCAL DENSITY OF STATES IN THE 
PRESENCE OF IMPURITY SCATTERING 

The effects of impurity scattering on the graphene 
LDOS have been extensively studiecP^in the past. It has 
been shown^I^'t^ that this gives rise to long- wavelength 
oscillations that decay as instead of the 1/r law ex- 
pected for conventional two-dimensional materials. Here 
we investigate how the merging of the two Dirac cones 
changes the form of these long- wavelength oscillations. 
We start this section with some zero-energy wavefunc- 
tion arguments along the lines of Ref. 14 that allow us to 



characterize the impurity state. Furthermore we perform 
a more detailed study using analytical (T-matrix approx- 
imation) and numerical (tight-binding) techniques. 

A. Wavefunction considerations 

1. The suhlattice symmetry 

Graphene honeycomb lattice contains two atoms per 
unit cell {A and 5), which allows one to define two 
sublattices. Moreover the Hamiltonian ([T]) only takes 
into account nearest neighbor hopping processes, and ne- 
glects hopping between sites belonging to the same sub- 
lattice, resulting in a bipartite system. For such systems 
a generic Hamiltonian takes the form: 

= T^t ) 

Without loss of generality, T is a Na x Nb block (not 
necessarily a square matrix), where N^^b) is the number 
of atoms in the A{B) sublattice, assuming there is only 
one electron per atom. Here we restrict ourself to Nb > 
A^^. Such a Hamiltonian anti-commutes with: 

I AT is the N X N identity matrix so that the unitary op- 
erator S always squares to +1, which defines a chiral 
symmetry: the sublattice symmetry. 

This fundamental symmetry implies a particle-hole 
symmetric spectrum, and includes the possibility of ex- 
istence of zero energy states, which transform into them- 
selves under the transformation S. As a consequence, 
they have null components on one sublattice. 

Moreover, as pointed out in Ref. IT9ti22| every finite 
bipartite lattice has an extra number of Nb — Na zero- 
energy eigenstates living on the sublattice 5, regard- 
less of the components of the block T. This is because 
the non-zero energy eigenstates appear in pairs: and 
tSIV^), and in order to form non-zero energy states, it is 
necessary to pair a localized state living on the sublattice 
A with another one living on B. As a number of Nb—Na 
zero modes living on the sublattice B are unable to sat- 
isfy this condition, they are stuck at zero energy, cannot 
hybridize to A states, and remain localized purely on B. 

2. Zero- energy impurity wavefunction 

In the presence of a single vacancy, Nb — Na = 1, and 
we have a single zero-mode impurity-state wavefunction. 
Here the fundamental point is that varying the param- 
eter t' does not change the structure of the matrix ([s]). 
Then the sublattice symmetry ensures that such a zero- 
mode does exist, even in the gapped phase {t' > 2). As 
a consequence, this zero-energy state is a good candidate 
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FIG. 3. The zero-energy wavefunction components for t' = at 
with a > 2. The wavefunction is zero for all sites for which no 
value is specified. The black square denotes the vacancy. The 
direction of anisotropy is the ^-direction. The two dashed 
lines are the boundaries of the upper and respectively lower 
half-planes. 



to characterize the Dirac cones merging in real space. 
In this section we study the form of its wavefunction, 
using simple arguments along the lines of Ref. IT, We 
already know that such a wavefunction has null com- 
ponents on the A sites, represented by the black disks 
in FigjSj and we need to determine its value on the B 
sublattice. In Ref. 14, the authors have determined the 
exact analytic form of the impurity wavefunction for an 
isotropic honeycomb lattice with a single vacancy. Their 
method consists in an appropriate matching of the zero 
modes of two semi-infinite and complementary graphene 
sheets. This is the method we generalize in what follows 
for anisotropic graphene. 

In Fig. |3] the two semi-infinite graphene sheets are 
defined such that their edges are orthogonal to the 
anisotropic direction, along which t' = at with a > 2. 
Here we have introduced an anisotropy parameter a, 
{a = 2 exactly at the merging) that allows us to explore 
the gapped phase beyond the merging point. The upper 
half-plane has a 'bearded' edge (as indicated by the up- 
per dashed line in Fig. [3|, whereas the lower half-plane 
has a zigzag edge. 

Let us first consider the lower half-plane terminated 
by the zigzag edge. The form of the edge states for a 
semi-infinite zigzag ribbon is well knowiP^^^for isotropic 
graphene. In a manner similar to that of Ref.'^S', the edge 
states for anisotropic graphene can be determined by im- 
posing the condition: |2cos(/c/2)| < a, where k is the 
momentum along the edge. While for isotropic graphene 
(a = 1) this condition is verified for 27r/3 < k < 47r/3, 
above the merging point {a > 2) such a condition is 
satisfied for all values of /c, < /c < 27r. Next, regard- 
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ing the complementary semi-infinite bearded plane, the 
condition becomes: |2cos(/c/2)| > a, which cannot be 
satisfied for any k when a > 2. The case a = 2 leads to 
/c = 0, associated to an extended state, and there are no 
allowed edge states in this limit. 

The condition that the impurity wavefunctions on the 
two semi-infinite planes match at the interface can be 
written as: 



p(r.E) for E=Q 



m,0 



^m+1,0 







(7) 



where bm,n {bin!n) corresponds to a given site of the lower 
(upper) half plane characterized by frn,n = Tn{a2 — ^i) — 
nai. The origin is defined to be on the B atom right 
below the vacancy in Fig. |3] The above relation is valid 
everywhere on the edges except for m = 0. Introducing 



^m,o = Z)/c ^^,06*^"^, the condition 
in terms of momentum as: 



can be rewritten 



a 



(0 ikm 



ik' \ Ak' m 



)e 



= 0. 



(8) 



A possible solution for the boundary conditions is b^l\ = 

1 with {) < k < 2tt and ^^.Tq = 0. As for the case of 
isotropic graphene studied in Ref. |14[ this corresponds to 
the edge solutions for two isolated complementary semi- 
infinite planes. Considering the lattice as infinite, the 
discrete sum in (|8| turns into an integral, and the impu- 
rity wavefunction can be written as: 



27r 



dk{-2/aY cos"(fc/2)e*'=(™+t ) 



(-1)^ 



(9) 



This approximation is valid for large distances when 
defining x = ao^/3{2m + ti)/2 and y = — n3ao/2. Most 
useful to compare with the results of the subsequent sec- 
tions is the behavior of the wavefunction along the direc- 
tion X = 0. Along this direction, the zero energy impu- 
rity state exhibits an exponential decay with the distance 
from the impurity in the gapped phase (a > 2), whereas 
it decays as l/^/y at the merging {a = 2). Moreover, 
if one evaluates the impurity wavefunction in the semi- 
metallic phase for 1 < a < 2, one can check that it still 
decreases as 1/r in both directions, albeit exhibiting a 
strong asymmetry between x and y. Hence the decay law 
of the zero-energy impurity states provides a real-space 
signature of the Dirac cone merging. 

Furthermore we can evaluate the amplitude of this im- 
purity state by hand, see Fig. [3) by searching for a de- 
caying wavefunction with null components everywhere in 
the semi-infinite bearded ribbon. The condition ^ be- 
comes Oib^^Q + + = and must be satisfied at each 
A site between the two dashed lines, except for the im- 
purity site. So the wavefunction has zero components 
along the zig-zag edge, except at the site situated right 
under the impurity, for which we take 6q q = 1. Then the 
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X (in units o f a) 



FIG. 4. Snapshot of the LDOS obtained using exact diagonal- 
ization. We plot the zero-energy impurity state slightly above 
the merging point, (a = 2.1) when a small gap opens in the 
spectrum. The highest-intensity site (in blue) corresponds to 
the B site right below the impurity. 



Hamiltonian rt5l implies that ab, 



m,l 



for all values of m. This leads to b^_^\ ^ = b^^ = —1/a 

and =0 for all other sites with n = 1. If we ex- 
tend this analysis to the subsequent rows, we obtain the 
impurity wavefunction values shown in Fig. [3] So above 
the merging point, this peculiar localized state describes 
electrons that are localized only in the lower-half plane 
of the graphene sheet with a single impurity. 
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B. Exact diagonalization 

In order to obtain the local density of states on the 
lattice in the presence of disorder one can diagonalize ex- 
actly or using numerical approximation the lattice tight- 
binding Hamiltonian. Here, since the systems we con- 
sider are not too large (around 1800 atoms) we base our- 
selves on an exact diagonalization technique. The lattice 
Hamiltonian is defined by: 



'^ = -E^^^-l^)01 + ^o|o)(o| 



(10) 



where \i) stands for the 2pz non- hybridized orbital cen- 
tered on site i. The impurity that we consider is a 
vacancy, which can be modeled by removing the corre- 
sponding atom from the lattice, or by taking an infinite 
value for Vq. We denote by \k) the eigenstate correspond- 
ing to the eigenvalue Ek. The eigenfunctions of (10) can 
be written as a linear combination of individual orbit als: 



\k) = Y.ck.\i), (11) 

i 

Cki = (12) 
The LDOS, corresponding to the number of available 
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states on a site i at energy E is then given by: 



PiiE) = J2\ckiffkiE) 



(13) 



where fk{E) = S{E — E^) is the Dirac delta function cen- 
tered on the eigenenergy E^ . While in an infinite system 
this procedure automatically yields a continuous energy 
spectrum, in a finite sample the spectrum is smoothed 
by taking //. to be a Gaussian or a Lorentzian. 

In Fig. |4] we show the LDOS obtained using this 
method at zero energy in the gapped phase {a = 2.1). 
This is in agreement with the zero-energy wavefunction 
described previously and depicted in Fig. |3] The result 
for the spatial dependence of the LDOS at a finite en- 
ergy is presented in Fig. [S] Note the strong asymmetry 
of the LDOS between the positive and negative values of 
y close to the impurity. While some of these features are 
consistent with the previous observations concerning the 
impurity-state wavefunctions, we also investigate them 
in more details (for example in what concerns their en- 
ergy dependence) in the next section, via the T-matrix 
approximation technique. 



C. T-matrix approximation 

The T-matrix approximation consists in a perturbative 
expansion of the Green's function to all orders in the 
impurity scattering, as shown in Fig. [5] In this paper we 
consider the case of a localized impurity V{f) = Vo5{r) 
situated on a sublattice A, for which V{q) is independent 
of momentum. We consider as impurity a vacancy, for 
which Vb becomes infinite. 

The expansion of the T-matrix in Fig. [5] is a geomet- 
ric series and the infinite summation of diagrams can be 
performed exactly: 



T{i^n) = [h-V 



I 

Jbz 



Sbz 



GS,i^n)]-^V (14) 



and V 



where Sbz is the area of the BZ, iuOn a Matsubara 
frequency, Go{k^iuJn) = [ioJnh — ^k\~^ is the unper- 
turbed Green's function, I2 the 2x2 identity matrix, 
>o 


We define Ap as the correction to the LDOS due to 
the impurity. According to Fig. [5j we have: 

AG{Ri, Rj, E) = G{Ri, Rj, E) - G^Ri - Rj, E) 

= G\Ri, E)T{E)G\-Rj, E) (15) 

The correspondence between the components of AG in 
the continuum and on the lattice is the following: 



AGo.p{ri,r2,E) = 



R^)AGa^{Rj,R^,E), (16) 



V 



V 



V 



FIG. 5. Diagrammatic perturbative expansion of the general- 
ized Green's function to all orders in the impurity potential. 



where (/>a(//3) is a carbon orbital, the a and P indices 
denote the sublattice, whereas i and j label the unit cell. 
The impurity correction to the LDOS is given by: 

Ap{r,E) = - ^Im[Tr{AG{r,r,E)}] (17) 



which yields in the momentum space: 
Ap{q,E) = 



(18) 



Jbz^bz 



Tr{AG{k + q, k, E) - AG%k, k^q,E)} 



1. Momentum dependence of the Fourier transform of the 
LDOS 

We focus first on the evaluation of the momentum de- 
pendence of Ap, corresponding to the measured Fourier 
transform of the LDOS. In Fig. [6] we plot this momentum 
dependence for uj = O.lbt and u = OM. The first col- 
umn corresponds to isotropic graphene {t' = t). As noted 
previousl};^^ the central circle (in red) corresponds to in- 
tranodal scattering, whereas the outer regions around the 
corners of the BZ correspond to internodal scattering. 
In the second column we consider an intermediate value 
of t' = 1.5t, while in the third column we consider the 
Dirac-cone merging limit f = 2t. We note that the outer 
regions disappear at the merging point for which intern- 
odal quasiparticle scattering no longer exists. Moreover, 
we note that the features corresponding to intranodal 
scattering, centered on the sites of the reciprocal lattice 
are strongly anisotropic, corresponding to the low-energy 
anisotropic semi-Dirac spectrum. 



2. Friedel oscillations in real space 

In what follows we focus on the real space form of the 
Friedel oscillations. While they can be evaluated numeri- 
cally for arbitrary energy and position using the formulae 
presented in the previous section, we can also obtain an 
analytical form of these oscillations in certain limits by 
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^Paa(Q:^) ^^Pbb(Q:^) ^Paa(q^^) ^^PebCQ:^) ^PaaCq^^) ^^PbeCq^^) 




-3 -2 -1 1 2 3 -3 -2 -1 1 2 3 -3 -2 -1 1 2 3 



{in units of a ^) {in units of a ^) q^ {in units of a ^) 

FIG. 6. Fourier transform of the LDOS correction Ap when t' — t (first column), t' — 1.5t (second column) and t' — 2t (third 
column). The energy is uj — 0.1 5t for the first row and uj — O.SOt for the second row. 



^PAA{'>^.^)+^PBB{'r.^) ^PAA{r.^)+^pBB{r,^) Ap^^(r,a;) +Ap55(f,a;) 




-20 -10 10 20 ■ -20 -10 10 20 -20 -10 10 20 



X {in units of a) x {in units of a) x {in units of a) 

FIG. 7. Correction of the LDOS Ap (in arbitrary units) for t' — t (first column), t' — 2t (second and third columns), obtained 
using the full T-matrix approximation (first and second columns) and the low-energy expansion (third column). The energy is 
u 0.15t. 
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performing an expansion of the Hamiltonian at low en- 
ergy, where the physics is dominated by the semi-Dirac 
spectrum around the M point. Using the expansion (|3| 
with A = 0, the unperturbed Green's function for x = 
can be rewritten as: 



G\0,y,Lo) 



_ f G%{0,y,u) G%{0,y,io) 



(19) 



with, 



G%(0,y,c.) = -i2-5/V3/2r(l/4)u;-V4yi/4j^(i)^^(^y) 
G%B(0,y,uj) = GAA{0,y,u;) 

G%{0,y,co) = -a2-3/4^3/2p(3/4)^i/4y-i/4^(i) (^y) 

TiAi2-5/V/2r(l/4).w3/4yl/4^(l)(^y) 

G%Ao,y,^) = -iA2-3/4^3/2r(3/4)^i/4y-i/4^(i)(^y) 



ai2-5/V3/2r(l/4)w3/4yl/4^W(^y) 



(20) 



where hI^^ are Hankel functions of the first kind, T is the 
Euler gamma function and A the conjugate of an arbi- 
trary phase factor A that depends on the basis we choose 
to write f{k); the value of any observable physical quan- 
tity should be independent of this choice^^. Note that on 
the right-hand-side of the above formulae we have chosen 
to denote the absolute value \y\ simply by y. Moreover, 
the =F signs correspond to a positive and respectively neg- 
ative value for y. The antisymmetric form of Gab and 
Gba is responsible for an anisotropy of the impurity- 
induced corrections to the LDOS on the B sublattice, as 
it will be described in more detail in what follows. Ac- 



cording to Eq. (17), the LDOS correction on each sub- 



lattice is given by: 

ApAA{0,y,uj) = -^Im[G%{0,y,uj)t{uj)G%{0,-y,uj)] 

ApBB{0,y,uj) = -^Im[G%A{0.y.^)tHG%{0,-y,uj)] 

(21) 

where t{uj) is the only non-nul component of the T- 
matrix. In the case of an infinite impurity potential (va- 
cancy), it takes the form: 



-in/A: 



,1/2 



(22) 



At this point, we can check that at zero energy our T- 
matrix calculations recover the same expression for the 
LDOS as the one obtained from the zero-energy wave- 
function considerations. In the limit a; ^ 0: 



G%{0,y,oj) ^ oj^/^ 
G%{0,y,uj) ^ A^9{y) 



(23) 



where is the Heaviside step function. The LDOS cor- 
rection then vanishes on the sublattice A. We stress that 
the sublattice symmetry implies that at zero energy the 
LDOS on the sublattice A is zero, whereas it behaves in 
the following manner on the sublattice B: 



ApBB{0,y,uj) 



oj-y) 



(24) 



This result is in agreement with the analysis of the zero- 
energy wavefunction. Remember that the impurity wave- 
function decays as \l ^Jy with the distance from the im- 
purity (cf. ([9|), which then leads to a Xjy decay for the 
LDOS. 

Now we turn back to the FO and evaluate the correc- 
tions to the LDOS using the corresponding expressions 
for the Green's function components in Eq. ([20|. The 
results are presented in Fig. [Sj We compare these results 
to a full evaluation of the T-matrix (without making the 
low-energy expansion), as well as with results obtained 
using the tight-binding method. As it can be seen in 
Fig. [8] all methods yield very similar results, which con- 
firms their accuracy for this type of calculation. We also 
note that the LDOS correction is asymmetric between 
the positive and negative values of y on the B sublattice, 
whereas it is symmetric on the A sublattice. 

To obtain the asymptotic expansion of the Friedel os- 
cillations we expand the Hankel functions for large values 
oiujy and we get: 



Apaa(0, y, ijJ) - cos ( h TT 1 



Apbb(0,^,cj) 



1/2 



cos 



y 



(25) 



The resulting FO decay as 1/y^ at large distances on 
both sublattices, slower than the typical inverse decay 
for a regular two-dimensional system, however their pe- 
riod is still proportional to l/cj. When summing the 
contribution of the two sublattices, the terms in 
which are dephased by a factor of tt, vanish. The FO 
are then described by the next leading correction which 
is non-zero only on the B sublattice: 



Ap(0,y,a;) 



,1/2 



y 



/2ujy 



(26) 



Here the minus/plus signs correspond to positive and re- 
spectively negative values of y. The long wavelength os- 
cillations thus decay following the usual 1/y law, different 
from the 1/r^ law corresponding to the intra- nodal scat- 
tering in typical graphene. Thus the transition from the 
1/r^ decay to a 1/r decay in the low-energy FO provides 
a real-space signature of the Dirac points merging. 

The Friedel oscillations along the perpendicular direc- 
tion {y = 0) cannot be evaluated analytically, however in 
the third column of Fig.|8]we plot the dependence of the 
Friedel oscillations as a function of x for y = 0. Note that 
the amplitude of the oscillations is greatly reduced with 
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FIG. 8. The LDOS correction as a function of position in the vicinity of the impurity. The second hne presents a series of zooms- 
in of the plots outUned on the first hne. In the first column we compare Ap obtained using the full T-matrix approximation 
to the one obtained by the tight-binding method for an energy cj = — 0.20t. Note that, consistent to the low energy expansion, 
the FO obey an inverse square-root decay for large values of y and are dephased by tt between the two sublattices. The second 
column presents a comparison between the correction to the LDOS Ap along the x — direction obtained by the full T-matrix 
approximation (full lines) and by the low energy expansion (dotted lines) for u = — 0.20t. In blue we plot the LDOS on the 
A sublattice comprising the impurity {y = 0) whereas in red the LDOS on the B sublattice. In the third column we plot the 
LDOS along the y — direction obtained by the full T-matrix approximation for uj — — 0.20t. The blue curve is obtained using 
the full T-matrix approximation while the black one is obtained in the continuum approximation. 
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respect of the oscillations in the y direction, consistent 
with the asymmetric shape of the impurity-state cloud, 
elongated in the ^/-direction. 

In Fig.[7|we also present a two-dimensional plot of the 
LDOS at a finite energy, obtained both by using the full 
T-matrix form and the low-energy expansion. Note that 
the behavior is very similar to that obtained using the 
tight-binding method described in the previous section. 

IV. CONCLUSION 

We have studied the LDOS in the presence of a simple 
impurity for an anisotropic graphene system at the Dirac- 
cone merging point. We have found that near this par- 
ticular point, the zero-energy impurity wavefunction and 
the Friedel oscillations in the LDOS exhibit very pecu- 
liar features. In particular, the decay length of the Friedel 
oscillations along the anisotropy direction and along the 
direction perpendicular to this direction are very differ- 
ent, yielding a very asymmetric impurity state in real 
space. The spatial dependence of the impurity state 
wavefunction allows us to clearly distinguish the semi- 
metallic phase (power-law decay of the wavefunction with 



the distance from the impurity) from the gapped phase 
(exponential decay). On the other hand, the semi-Dirac 
spectrum near the merging point induces a change of 
the decay laws in the Friedel oscillations from a inverse- 
square law (l/r^) below the transition to an inverse-linear 
law (1/r) exactly at the transition. At low energy this 
provides a real-space signature of the topological transi- 
tion. 

The agreement between the methods that we have 
used, the analytical T-matrix approximation, the numeri- 
cal tight-binding exact diagonalization and the wavefunc- 
tion considerations, is remarkable, proving the accuracy 
of these methods to describe the LDOS in the presence 
of disorder in a generic two-dimensional system. 
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